Community Analysis
Daniyal Tariq
2024-04-30
##Goals
Load in phyloseq data with rooted tree. Evaluate sequencing depth and remove sample. Normalize the read counts between samples. Calculate community dissimilarities. Numbers between 0 and 1. If 0, completely similar versus if they are 1, then they’re completely dissimilar. Sorensen: Shared Species as a binary value: Abundance-unweighted Bray-Curtis: Shared Abundant species: Abundance-weighted (Abundance-)Weighted UNIFRAC: Consider Abundant Species and where they fall on the tree Visualize the community data with two unconstrained Ordinations: PCoA: Linear Method. Eigenvalue = how much variation is explained by each axis. Choose to view axis 1, 2, 3, etc. and plot them together. NMDS: Non-linear. Smush multiple Dimensions into 2 or 3 axes. Need to report Stress value (ideally <0.15). Run statistics with PERMANOVA and betadispR.
##Set seed
set.seed(2443)##Load libraries
#install.packages("vegan")
pacman::p_load(tidyverse, devtools, phyloseq, patchwork, vegan,
install = FALSE)
library(gridExtra)##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combine# Load Lakestie colors
lakesite_colors <- c("110" = "dodgerblue4","M15" = "dodgerblue2","M45" = "#D9CC3C","MLB" = "#A0E0BA")##Load Data
# Load in rooted phylogenetic tree!
load("data/03_Phylogenetic_Tree/phytree_preprocessed_physeq.RData")##Explore Read Counts
# Calculate the total number of reads per sample.
raw_TotalSeqs_df_combined <-
midroot_physeq %>%
# calculate the sample read sums
sample_sums() %>%
data.frame()
# name the column
colnames(raw_TotalSeqs_df_combined)[1] <- "TotalSeqs"
raw_TotalSeqs_df_DNA <-
midroot_physeq_DNA %>%
# calculate the sample read sums
sample_sums() %>%
data.frame()
# name the column
colnames(raw_TotalSeqs_df_DNA)[1] <- "TotalSeqs"
raw_TotalSeqs_df_RNA <-
midroot_physeq_RNA %>%
# calculate the sample read sums
sample_sums() %>%
data.frame()
# name the column
colnames(raw_TotalSeqs_df_RNA)[1] <- "TotalSeqs"
head(raw_TotalSeqs_df_combined)## TotalSeqs
## 110D2W115D 13235
## 110D2W115R 12953
## 110D2W415D 14765
## 110D2W415R 11587
## 110D2W515D 18867
## 110D2W515R 12432head(raw_TotalSeqs_df_DNA)## TotalSeqs
## 110D2W115D 13235
## 110D2W415D 14765
## 110D2W515D 18867
## 110D2W615D 9289
## 110D2W715D 14757
## 110D2W815D 36720head(raw_TotalSeqs_df_RNA)## TotalSeqs
## 110D2W115R 12953
## 110D2W415R 11587
## 110D2W515R 12432
## 110D2W615R 16887
## 110D2W715R 14711
## 110D2W815R 10574# make histogram of raw reads
raw_TotalSeqs_df_combined %>%
ggplot(aes(x = TotalSeqs)) +
geom_histogram(bins = 50) +
scale_x_continuous(limits = c(0, 10000)) +
labs(title = "Raw Combined Sequencing Depth Distribution") +
theme_classic()## Warning: Removed 74 rows containing non-finite outside the scale range
## (`stat_bin()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_bar()`).
raw_TotalSeqs_df_DNA %>%
ggplot(aes(x = TotalSeqs)) +
geom_histogram(bins = 50) +
scale_x_continuous(limits = c(0, 10000)) +
labs(title = "Raw DNA Sequencing Depth Distribution") +
theme_classic()## Warning: Removed 43 rows containing non-finite outside the scale range (`stat_bin()`).
## Removed 2 rows containing missing values or values outside the scale range
## (`geom_bar()`).
raw_TotalSeqs_df_RNA %>%
ggplot(aes(x = TotalSeqs)) +
geom_histogram(bins = 50) +
scale_x_continuous(limits = c(0, 10000)) +
labs(title = "Raw RNA Sequencing Depth Distribution") +
theme_classic()## Warning: Removed 31 rows containing non-finite outside the scale range (`stat_bin()`).
## Removed 2 rows containing missing values or values outside the scale range
## (`geom_bar()`).
#What are the minimum number of sequences for each rooted physeq object?
midroot_physeq %>% sample_sums() %>% min()## [1] 2142midroot_physeq_DNA %>% sample_sums() %>% min()## [1] 2582midroot_physeq_RNA %>% sample_sums() %>% min()## [1] 2142#The minimum number of sequences for DNA samples = 2582 and for RNA samples = 2142#In order to normailze the read counts we first need to import a scale read function from http://deneflab.github.io/MicrobeMiseq/
matround <- function(x){trunc(x+0.5)}
scale_reads <- function(physeq, n = min(sample_sums(physeq)), round = "round") {
# transform counts to n
physeq.scale <- transform_sample_counts(physeq, function(x) {(n * x/sum(x))})
# Pick the rounding functions
if (round == "floor"){
otu_table(physeq.scale) <- floor(otu_table(physeq.scale))
} else if (round == "round"){
otu_table(physeq.scale) <- round(otu_table(physeq.scale))
} else if (round == "matround"){
otu_table(physeq.scale) <- matround(otu_table(physeq.scale))
}
# Prune taxa and return new phyloseq object
physeq.scale <- prune_taxa(taxa_sums(physeq.scale) > 0, physeq.scale)
return(physeq.scale)
}#Scale the reads using the above function, not we are scaling each object by their minimum number for sequences.
# Scale reads combined DNA/RNA by the above function
scaled_rooted_physeq <-
midroot_physeq %>%
scale_reads(round = "matround")
# Scale reads DNA samples by the above function
scaled_rooted_physeq_DNA <-
midroot_physeq_DNA %>%
scale_reads(round = "matround")
# Scale reads RNA samples by the above function
scaled_rooted_physeq_RNA <-
midroot_physeq_RNA %>%
scale_reads(round = "matround")
# Calculate combined DNA//RNA read depth
scaled_TotalSeqs_df <-
scaled_rooted_physeq %>%
sample_sums() %>%
data.frame()
colnames(scaled_TotalSeqs_df)[1] <-"TotalSeqs"
head(scaled_TotalSeqs_df) #Inspect to see if the scaling occured## TotalSeqs
## 110D2W115D 2143
## 110D2W115R 2130
## 110D2W415D 2137
## 110D2W415R 2140
## 110D2W515D 2146
## 110D2W515R 2147# Calculate DNA samples read depth
scaled_TotalSeqs_DNA_df <-
scaled_rooted_physeq_DNA %>%
sample_sums() %>%
data.frame()
colnames(scaled_TotalSeqs_DNA_df)[1] <-"TotalSeqs"
head(scaled_TotalSeqs_DNA_df) #Inspect to see if the scaling occured## TotalSeqs
## 110D2W115D 2581
## 110D2W415D 2583
## 110D2W515D 2584
## 110D2W615D 2586
## 110D2W715D 2583
## 110D2W815D 2568# Calculate RNA samples read depth
scaled_TotalSeqs_RNA_df <-
scaled_rooted_physeq_RNA %>%
sample_sums() %>%
data.frame()
colnames(scaled_TotalSeqs_RNA_df)[1] <-"TotalSeqs"
# Inspect
head(scaled_TotalSeqs_RNA_df)## TotalSeqs
## 110D2W115R 2130
## 110D2W415R 2140
## 110D2W515R 2147
## 110D2W615R 2138
## 110D2W715R 2142
## 110D2W815R 2135#Check the range of the data
#For Combined,DNA and RNA
min_seqs_combined <- min(scaled_TotalSeqs_df$TotalSeqs);min_seqs_combined## [1] 2117min_seqs_DNA <- min(scaled_TotalSeqs_DNA_df$TotalSeqs); min_seqs_DNA## [1] 2568min_seqs_RNA <- min(scaled_TotalSeqs_RNA_df$TotalSeqs); min_seqs_RNA## [1] 2129max_seqs_combined <- max(scaled_TotalSeqs_df$TotalSeqs);max_seqs_combined## [1] 2173max_seqs_DNA <- max(scaled_TotalSeqs_DNA_df$TotalSeqs); max_seqs_DNA## [1] 2595max_seqs_RNA <- max(scaled_TotalSeqs_RNA_df$TotalSeqs); max_seqs_RNA## [1] 2173#Range
max_seqs_combined - min_seqs_combined## [1] 56max_seqs_DNA - min_seqs_DNA## [1] 27max_seqs_RNA - min_seqs_RNA## [1] 44# Plot Histogram
scaled_TotalSeqs_df %>%
ggplot(aes(x = TotalSeqs)) +
geom_histogram(bins = 50) +
scale_x_continuous(limits = c(0, 10000)) +
labs(title = "Scaled Sequencing Depth at 2142") +
theme_classic()## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_bar()`).
scaled_TotalSeqs_DNA_df %>%
ggplot(aes(x = TotalSeqs)) +
geom_histogram(bins = 50) +
scale_x_continuous(limits = c(0, 10000)) +
labs(title = "Scaled Sequencing Depth of DNA only samples at 2582") +
theme_classic()## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_bar()`).
scaled_TotalSeqs_RNA_df %>%
ggplot(aes(x = TotalSeqs)) +
geom_histogram(bins = 50) +
scale_x_continuous(limits = c(0, 10000)) +
labs(title = "Scaled Sequencing Depth of RNA only samples at 2142") +
theme_classic()## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_bar()`).
#We can now start running PCoA with respect to Soreson and Bray-Curtis
# Calculate sorensen dissimilarity: Abundance-unweighted of shared taxa
# Running PCoA with limnion as the shape
scaled_soren_pcoa_DNA <-
ordinate(
physeq = scaled_rooted_physeq_DNA,
method = "PCoA",
distance = "bray", binary = TRUE)
#str(scaled_soren_pcoa_DNA)
# Plot the ordination
soren_limnion_pcoa_DNA <- plot_ordination(
physeq = scaled_rooted_physeq_DNA,
ordination = scaled_soren_pcoa_DNA,
color = "lakesite",
shape = "limnion",
title = "Sorensen PCoA All Taxa (DNA)") +
geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
scale_shape_manual(values = c(15,16,17,18, 19, 20)) +
scale_color_manual(values = lakesite_colors) +
theme_bw()
# Calculate sorensen dissimilarity: Abundance-unweighted of shared taxa
scaled_soren_pcoa_RNA <-
ordinate(
physeq = scaled_rooted_physeq_RNA,
method = "PCoA",
distance = "bray", binary = TRUE)
#str(scaled_soren_pcoa_DNA)
# Plot the ordination
soren_limnion_pcoa_RNA <- plot_ordination(
physeq = scaled_rooted_physeq_RNA,
ordination = scaled_soren_pcoa_RNA,
color = "lakesite",
shape = "limnion",
title = "Sorensen PCoA Active Taxa (RNA)") +
geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
scale_shape_manual(values = c(15,16,17,18, 19, 20)) +
scale_color_manual(values = lakesite_colors) +
theme_bw()
# Show the plot
grid.arrange(soren_limnion_pcoa_DNA,soren_limnion_pcoa_RNA, ncol=2) 
#Running PCoA with season as the shape
# Plot the ordination
soren_season_pcoa_DNA <- plot_ordination(
physeq = scaled_rooted_physeq_DNA,
ordination = scaled_soren_pcoa_DNA,
color = "lakesite",
shape = "season",
title = "Sorensen PCoA DNA") +
geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
scale_shape_manual(values = c(15,16,17,18, 19, 20)) +
scale_color_manual(values = lakesite_colors) +
theme_bw()
# Plot the ordination
soren_season_pcoa_RNA <- plot_ordination(
physeq = scaled_rooted_physeq_RNA,
ordination = scaled_soren_pcoa_RNA,
color = "lakesite",
shape = "season",
title = "Sorensen PCoA RNA") +
geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
scale_shape_manual(values = c(15,16,17,18, 19, 20)) +
scale_color_manual(values = lakesite_colors) +
theme_bw()
# Show the plot
grid.arrange(soren_season_pcoa_DNA,soren_season_pcoa_RNA,ncol=2)
The PCoAs explain a fair amount of variation for RNA/DNA about 47% for DNA and 40.3% for RNA. The clustering seems quite similar between the DNA and RNA samples this indicates with the MLB(surface) samples seem to cluster in two regions. The liminions do not seem to match indicating that the the variation may not be due to the liminion. The other samples seem to be spread quite wide with no clear pattern arising. Though the DNA and RNA seem to cluster similarly a point of note is that there seems to be a larger difference in the RNA surface MLB cluster on the second axis indicating that there is a difference in active taxa relative to total taxa, this does not seem to be limnion dependant.
# Calculate the BC distance
scaled_BC_pcoa_DNA <-
ordinate(
physeq = scaled_rooted_physeq_DNA,
method = "PCoA",
distance = "bray")
scaled_BC_pcoa_RNA <-
ordinate(
physeq = scaled_rooted_physeq_RNA,
method = "PCoA",
distance = "bray")
# Plot the PCoA
bray_pcoa_limnion_DNA <-
plot_ordination(
physeq = scaled_rooted_physeq_DNA,
ordination = scaled_BC_pcoa_DNA,
color = "lakesite",
shape = "limnion",
title = "Bray-Curtis PCoA DNA Limnion") +
geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
scale_shape_manual(values = c(15,16,17,18, 19, 20,21)) +
scale_color_manual(values = c(lakesite_colors)) +
theme_bw()
bray_pcoa_limnion_RNA <-
plot_ordination(
physeq = scaled_rooted_physeq_RNA,
ordination = scaled_BC_pcoa_RNA,
color = "lakesite",
shape = "limnion",
title = "Bray-Curtis PCoA RNA Limnion") +
geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
scale_shape_manual(values = c(15,16,17,18, 19, 20,21)) +
scale_color_manual(values = c(lakesite_colors)) +
theme_bw()
bray_pcoa_season_DNA <-
plot_ordination(
physeq = scaled_rooted_physeq_DNA,
ordination = scaled_BC_pcoa_DNA,
color = "lakesite",
shape= "season",
title = "Bray-Curtis PCoA DNA Season") +
geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
scale_shape_manual(values = c(15,16,17,18, 19, 20,21)) +
scale_color_manual(values = c(lakesite_colors)) +
theme_bw()
bray_pcoa_season_RNA <-
plot_ordination(
physeq = scaled_rooted_physeq_RNA,
ordination = scaled_BC_pcoa_RNA,
color = "lakesite",
shape= "season",
title = "Bray-Curtis PCoA RNA Season") +
geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
scale_shape_manual(values = c(15,16,17,18, 19, 20,21)) +
scale_color_manual(values = c(lakesite_colors)) +
theme_bw()
#Show limnion plot
grid.arrange(bray_pcoa_limnion_DNA,
bray_pcoa_limnion_RNA, ncol=2)
#Show season plot
grid.arrange(bray_pcoa_season_DNA,
bray_pcoa_season_RNA, ncol=2) 
The PCOAs explain a high percentage of the variance for both the DNA (~45%) and RNA(~40%) samples. The top plots look at the limnion dependence of the abundant taxa and we can see that the points aren’t clustering based on shape. However, we see that the points cluster closely based on both shape (representing season) and color (representing lakesite) for the season plots. We see that the surface taxa once again cluster distinctly along the first axis for all these plots, indicating that they are quite different from the other depths that were sampled.
Weighted Unifrac
scaled_wUNI_pcoa_combined <-
ordinate(
physeq = scaled_rooted_physeq,
method = "PCoA",
distance = "wunifrac")
scaled_wUNI_pcoa_DNA <-
ordinate(
physeq = scaled_rooted_physeq_DNA,
method = "PCoA",
distance = "wunifrac")
scaled_wUNI_pcoa_RNA <-
ordinate(
physeq = scaled_rooted_physeq_RNA,
method = "PCoA",
distance = "wunifrac")
# Plot the PCoA
wUNI_station_pcoa_combined_season <-
plot_ordination(
physeq = scaled_rooted_physeq,
ordination = scaled_wUNI_pcoa_combined,
color = "lakesite",
shape = "season",
title = "Weighted Unifrac PCoA Combined season") +
geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
scale_shape_manual(values = c(15,16,17,18, 19, 20)) +
scale_color_manual(values = c(lakesite_colors)) +
theme_bw()
wUNI_station_pcoa_combined_limnion <-
plot_ordination(
physeq = scaled_rooted_physeq,
ordination = scaled_wUNI_pcoa_combined,
color = "lakesite",
shape = "limnion",
title = "Weighted Unifrac PCoA Combined limnion") +
geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
scale_shape_manual(values = c(15,16,17,18, 19, 20)) +
scale_color_manual(values = c(lakesite_colors)) +
theme_bw()
wUNI_station_pcoa_DNA_season <-
plot_ordination(
physeq = scaled_rooted_physeq_DNA,
ordination = scaled_wUNI_pcoa_DNA,
color = "lakesite",
shape = "season",
title = "Weighted Unifrac PCoA DNA season") +
geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
scale_shape_manual(values = c(15,16,17,18, 19, 20)) +
scale_color_manual(values = c(lakesite_colors)) +
theme_bw()
wUNI_station_pcoa_DNA_limnion <-
plot_ordination(
physeq = scaled_rooted_physeq_DNA,
ordination = scaled_wUNI_pcoa_DNA,
color = "lakesite",
shape = "limnion",
title = "Weighted Unifrac PCoA DNA limnion") +
geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
scale_shape_manual(values = c(15,16,17,18, 19, 20)) +
scale_color_manual(values = c(lakesite_colors)) +
theme_bw()
wUNI_station_pcoa_RNA_season <-
plot_ordination(
physeq = scaled_rooted_physeq_RNA,
ordination = scaled_wUNI_pcoa_RNA,
color = "lakesite",
shape = "season",
title = "Weighted Unifrac PCoA RNA season") +
geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
scale_shape_manual(values = c(15,16,17,18, 19, 20)) +
scale_color_manual(values = c(lakesite_colors)) +
theme_bw()
wUNI_station_pcoa_RNA_limnion <-
plot_ordination(
physeq = scaled_rooted_physeq_RNA,
ordination = scaled_wUNI_pcoa_RNA,
color = "lakesite",
shape = "limnion",
title = "Weighted Unifrac PCoA RNA limnion") +
geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
scale_shape_manual(values = c(15,16,17,18, 19, 20)) +
scale_color_manual(values = c(lakesite_colors)) +
theme_bw()
grid.arrange(wUNI_station_pcoa_combined_season,wUNI_station_pcoa_combined_limnion,ncol=2)
#Difference between seasons
grid.arrange(wUNI_station_pcoa_DNA_season, wUNI_station_pcoa_RNA_season,ncol=2)
#Difference between depths
grid.arrange(wUNI_station_pcoa_DNA_limnion, wUNI_station_pcoa_RNA_limnion,ncol=2)
The PCOAs explain a high percentage of the variance for both the DNA (~45%) and RNA(~47%) samples. We do not see a limnion or season dependence for the weighted uniFRAC distance. We did note one interesting feature, however. For the 110 site, there seems to be a cluster in the DNA samples showing that the total taxa present at this site are quite evolutionary similar. The active taxa, however, which we are inferring from the RNA samples show no clustering indicating that they are phylogenetically different.
# Calculate the Weighted Unifrac distance
scaled_wUNI_nmds_combined <-
ordinate(
physeq = scaled_rooted_physeq,
method = "NMDS",
distance = "wunifrac")## Run 0 stress 0.1418455
## Run 1 stress 0.14491
## Run 2 stress 0.1615657
## Run 3 stress 0.1562387
## Run 4 stress 0.166838
## Run 5 stress 0.1619733
## Run 6 stress 0.1494274
## Run 7 stress 0.1477715
## Run 8 stress 0.1550695
## Run 9 stress 0.147955
## Run 10 stress 0.1542195
## Run 11 stress 0.1449249
## Run 12 stress 0.1542195
## Run 13 stress 0.1668388
## Run 14 stress 0.1449098
## Run 15 stress 0.1418455
## ... New best solution
## ... Procrustes: rmse 1.651949e-05 max resid 0.0001469256
## ... Similar to previous best
## Run 16 stress 0.1634479
## Run 17 stress 0.1604705
## Run 18 stress 0.147952
## Run 19 stress 0.1600244
## Run 20 stress 0.1418454
## ... New best solution
## ... Procrustes: rmse 3.563792e-05 max resid 0.0003311476
## ... Similar to previous best
## *** Best solution repeated 1 timesscaled_wUNI_nmds_DNA <-
ordinate(
physeq = scaled_rooted_physeq_DNA,
method = "NMDS",
distance = "wunifrac")## Run 0 stress 0.1925474
## Run 1 stress 0.1985646
## Run 2 stress 0.2141103
## Run 3 stress 0.1985623
## Run 4 stress 0.1925474
## ... Procrustes: rmse 4.622689e-06 max resid 1.70631e-05
## ... Similar to previous best
## Run 5 stress 0.2515052
## Run 6 stress 0.1925474
## ... Procrustes: rmse 1.281947e-05 max resid 8.138397e-05
## ... Similar to previous best
## Run 7 stress 0.1926273
## ... Procrustes: rmse 0.004161513 max resid 0.02407351
## Run 8 stress 0.1925474
## ... Procrustes: rmse 3.626435e-05 max resid 0.0002383358
## ... Similar to previous best
## Run 9 stress 0.2102099
## Run 10 stress 0.1926274
## ... Procrustes: rmse 0.004159529 max resid 0.02406611
## Run 11 stress 0.198565
## Run 12 stress 0.1985652
## Run 13 stress 0.2101843
## Run 14 stress 0.1926273
## ... Procrustes: rmse 0.004166649 max resid 0.02412437
## Run 15 stress 0.1925474
## ... Procrustes: rmse 1.676562e-05 max resid 0.0001045908
## ... Similar to previous best
## Run 16 stress 0.1985644
## Run 17 stress 0.1985654
## Run 18 stress 0.1925474
## ... Procrustes: rmse 6.392559e-06 max resid 2.872485e-05
## ... Similar to previous best
## Run 19 stress 0.2451172
## Run 20 stress 0.1925474
## ... Procrustes: rmse 7.540351e-06 max resid 4.720985e-05
## ... Similar to previous best
## *** Best solution repeated 6 timesscaled_wUNI_nmds_RNA <-
ordinate(
physeq = scaled_rooted_physeq_RNA,
method = "NMDS",
distance = "wunifrac")## Run 0 stress 0.1645868
## Run 1 stress 0.2660364
## Run 2 stress 0.1645868
## ... New best solution
## ... Procrustes: rmse 4.067466e-06 max resid 2.252153e-05
## ... Similar to previous best
## Run 3 stress 0.1645868
## ... Procrustes: rmse 4.083904e-06 max resid 1.905353e-05
## ... Similar to previous best
## Run 4 stress 0.1645868
## ... Procrustes: rmse 2.722264e-06 max resid 1.323778e-05
## ... Similar to previous best
## Run 5 stress 0.1645868
## ... Procrustes: rmse 1.98128e-06 max resid 1.034146e-05
## ... Similar to previous best
## Run 6 stress 0.1645868
## ... Procrustes: rmse 2.849117e-06 max resid 1.05118e-05
## ... Similar to previous best
## Run 7 stress 0.1645868
## ... Procrustes: rmse 5.445529e-06 max resid 3.531163e-05
## ... Similar to previous best
## Run 8 stress 0.1645868
## ... Procrustes: rmse 2.938561e-06 max resid 1.758468e-05
## ... Similar to previous best
## Run 9 stress 0.1645868
## ... Procrustes: rmse 5.152021e-06 max resid 2.938458e-05
## ... Similar to previous best
## Run 10 stress 0.1645868
## ... Procrustes: rmse 3.21662e-06 max resid 1.860854e-05
## ... Similar to previous best
## Run 11 stress 0.2444954
## Run 12 stress 0.1645868
## ... Procrustes: rmse 4.164196e-06 max resid 2.150385e-05
## ... Similar to previous best
## Run 13 stress 0.1645868
## ... New best solution
## ... Procrustes: rmse 1.422965e-06 max resid 4.641084e-06
## ... Similar to previous best
## Run 14 stress 0.1645868
## ... Procrustes: rmse 3.82778e-06 max resid 1.75597e-05
## ... Similar to previous best
## Run 15 stress 0.2249074
## Run 16 stress 0.1645868
## ... Procrustes: rmse 3.358032e-06 max resid 2.109151e-05
## ... Similar to previous best
## Run 17 stress 0.1645868
## ... Procrustes: rmse 5.317784e-06 max resid 2.250939e-05
## ... Similar to previous best
## Run 18 stress 0.1645868
## ... Procrustes: rmse 1.990947e-06 max resid 7.963673e-06
## ... Similar to previous best
## Run 19 stress 0.1645868
## ... Procrustes: rmse 7.68143e-06 max resid 4.757098e-05
## ... Similar to previous best
## Run 20 stress 0.1645868
## ... Procrustes: rmse 2.829496e-06 max resid 1.629448e-05
## ... Similar to previous best
## *** Best solution repeated 7 times# Plot the PCoA
wUNI_station_nmds_combined_season <-
plot_ordination(
physeq = scaled_rooted_physeq,
ordination = scaled_wUNI_nmds_combined,
color = "lakesite",
shape = "season",
title = "Weighted Unifrac NMDS combined season") +
geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
scale_shape_manual(values = c(15,16,17)) +
scale_color_manual(values = c(lakesite_colors)) +
theme_bw()
wUNI_station_nmds_combined_limnion <-
plot_ordination(
physeq = scaled_rooted_physeq,
ordination = scaled_wUNI_nmds_combined,
color = "lakesite",
shape = "limnion",
title = "Weighted Unifrac NMDS combined limnion") +
geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
scale_shape_manual(values = c(15,16,17)) +
scale_color_manual(values = c(lakesite_colors)) +
theme_bw()
wUNI_station_nmds_DNA_season <-
plot_ordination(
physeq = scaled_rooted_physeq_DNA,
ordination = scaled_wUNI_nmds_DNA,
color = "lakesite",
shape = "season",
title = "Weighted Unifrac NMDS DNA season") +
geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
scale_shape_manual(values = c(15,16,17)) +
scale_color_manual(values = c(lakesite_colors)) +
theme_bw()
wUNI_station_nmds_DNA_limnion <-
plot_ordination(
physeq = scaled_rooted_physeq_DNA,
ordination = scaled_wUNI_nmds_DNA,
color = "lakesite",
shape = "limnion",
title = "Weighted Unifrac NMDS DNA limnion") +
geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
scale_shape_manual(values = c(15,16,17)) +
scale_color_manual(values = c(lakesite_colors)) +
theme_bw()
wUNI_station_nmds_RNA_season <-
plot_ordination(
physeq = scaled_rooted_physeq_RNA,
ordination = scaled_wUNI_nmds_RNA,
color = "lakesite",
shape = "season",
title = "Weighted Unifrac NMDS RNA season") +
geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
scale_shape_manual(values = c(15,16,17)) +
scale_color_manual(values = c(lakesite_colors)) +
theme_bw()
wUNI_station_nmds_RNA_limnion <-
plot_ordination(
physeq = scaled_rooted_physeq_RNA,
ordination = scaled_wUNI_nmds_RNA,
color = "lakesite",
shape = "limnion",
title = "Weighted Unifrac NMDS RNA limnion") +
geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
scale_shape_manual(values = c(15,16,17)) +
scale_color_manual(values = c(lakesite_colors)) +
theme_bw()
grid.arrange(wUNI_station_nmds_combined_season, wUNI_station_nmds_combined_limnion, ncol=2)
grid.arrange(wUNI_station_nmds_DNA_season, wUNI_station_nmds_RNA_season,ncol=2)
grid.arrange(wUNI_station_nmds_DNA_limnion, wUNI_station_nmds_RNA_limnion,ncol=2)
These NMDS plots don’t reveal any salient features, other than clustering based on lakesite.
##Statistical Significance Testing
# Calculate all three of the distance matrices
scaled_sorensen_dist_combined <- phyloseq::distance(scaled_rooted_physeq, method = "bray", binary = TRUE)
scaled_sorensen_dist_DNA <- phyloseq::distance(scaled_rooted_physeq_DNA, method = "bray", binary = TRUE)
scaled_sorensen_dist_RNA <- phyloseq::distance(scaled_rooted_physeq_RNA, method = "bray", binary = TRUE)
scaled_bray_dist_combined <- phyloseq::distance(scaled_rooted_physeq, method = "bray")
scaled_bray_dist_DNA <- phyloseq::distance(scaled_rooted_physeq_DNA, method = "bray")
scaled_bray_dist_RNA <- phyloseq::distance(scaled_rooted_physeq_RNA, method = "bray")
scaled_wUnifrac_dist_combined <- phyloseq::distance(scaled_rooted_physeq, method = "wunifrac")
scaled_wUnifrac_dist_DNA <- phyloseq::distance(scaled_rooted_physeq_DNA, method = "wunifrac")
scaled_wUnifrac_dist_RNA <- phyloseq::distance(scaled_rooted_physeq_RNA, method = "wunifrac")
# make a data frame from the sample_data
# All distance matrices will be the same metadata because they
# originate from the same phyloseq object.
metadata_combined <- data.frame(sample_data(scaled_rooted_physeq))
metadata_DNA <- data.frame(sample_data(scaled_rooted_physeq_DNA))
metadata_RNA <- data.frame(sample_data(scaled_rooted_physeq_RNA))
# Adonis test
# In this example we are testing the hypothesis that the five stations
# that were collected have different centroids in the ordination space
# for each of the dissimilarity metrics, we are using a discrete variable
adonis2(scaled_sorensen_dist_combined ~ lakesite, data = metadata_combined)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = scaled_sorensen_dist_combined ~ lakesite, data = metadata_combined)
## Df SumOfSqs R2 F Pr(>F)
## lakesite 3 5.8697 0.23886 10.983 0.001 ***
## Residual 105 18.7043 0.76114
## Total 108 24.5740 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(scaled_bray_dist_combined ~ lakesite, data = metadata_combined)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = scaled_bray_dist_combined ~ lakesite, data = metadata_combined)
## Df SumOfSqs R2 F Pr(>F)
## lakesite 3 4.4486 0.16638 6.9855 0.001 ***
## Residual 105 22.2892 0.83362
## Total 108 26.7379 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(scaled_wUnifrac_dist_combined ~ lakesite, data = metadata_combined)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = scaled_wUnifrac_dist_combined ~ lakesite, data = metadata_combined)
## Df SumOfSqs R2 F Pr(>F)
## lakesite 3 0.5782 0.08178 3.1171 0.002 **
## Residual 105 6.4923 0.91822
## Total 108 7.0705 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(scaled_sorensen_dist_DNA ~ lakesite, data = metadata_DNA)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = scaled_sorensen_dist_DNA ~ lakesite, data = metadata_DNA)
## Df SumOfSqs R2 F Pr(>F)
## lakesite 3 3.2802 0.32063 8.1804 0.001 ***
## Residual 52 6.9503 0.67937
## Total 55 10.2305 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(scaled_bray_dist_DNA ~ lakesite, data = metadata_DNA)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = scaled_bray_dist_DNA ~ lakesite, data = metadata_DNA)
## Df SumOfSqs R2 F Pr(>F)
## lakesite 3 2.5944 0.30663 7.6653 0.001 ***
## Residual 52 5.8667 0.69337
## Total 55 8.4611 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(scaled_wUnifrac_dist_DNA ~ lakesite, data = metadata_DNA)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = scaled_wUnifrac_dist_DNA ~ lakesite, data = metadata_DNA)
## Df SumOfSqs R2 F Pr(>F)
## lakesite 3 0.28642 0.17106 3.5768 0.001 ***
## Residual 52 1.38801 0.82894
## Total 55 1.67444 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(scaled_sorensen_dist_RNA ~ lakesite, data = metadata_RNA)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = scaled_sorensen_dist_RNA ~ lakesite, data = metadata_RNA)
## Df SumOfSqs R2 F Pr(>F)
## lakesite 3 3.1514 0.25969 5.7294 0.001 ***
## Residual 49 8.9840 0.74031
## Total 52 12.1354 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(scaled_bray_dist_RNA ~ lakesite, data = metadata_RNA)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = scaled_bray_dist_RNA ~ lakesite, data = metadata_RNA)
## Df SumOfSqs R2 F Pr(>F)
## lakesite 3 2.7364 0.23186 4.9302 0.001 ***
## Residual 49 9.0656 0.76814
## Total 52 11.8020 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(scaled_wUnifrac_dist_RNA ~ lakesite, data = metadata_RNA)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = scaled_wUnifrac_dist_RNA ~ lakesite, data = metadata_RNA)
## Df SumOfSqs R2 F Pr(>F)
## lakesite 3 0.420 0.18049 3.5973 0.001 ***
## Residual 49 1.907 0.81951
## Total 52 2.327 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Using lakesite as our primary explanatory variable, we see that we see statistical significance for all our dissimilarity metrics when we conduct a PERMANOVA analysis. This indicates that the community composition is different at a statistically significant level between each of the different depths that were sampled and this is true for the total taxa (indicated by the DNA samples) as well as the active taxa (indicated by the RNA samples).
##BetaDispR
# Homogeneity of Disperson test with beta dispr
# Sorensen
beta_soren_station_combined <- betadisper(scaled_sorensen_dist_combined, metadata_combined$lakesite)
permutest(beta_soren_station_combined)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 0.14349 0.047831 5.6975 999 0.004 **
## Residuals 105 0.88148 0.008395
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Bray-curtis
beta_bray_station_combined <- betadisper(scaled_bray_dist_combined, metadata_combined$lakesite)
permutest(beta_bray_station_combined)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 0.04498 0.014993 1.7103 999 0.153
## Residuals 105 0.92043 0.008766# Weighted Unifrac
beta_unifrac_station_combined <- betadisper(scaled_wUnifrac_dist_combined, metadata_combined$lakesite)
permutest(beta_unifrac_station_combined)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 0.001152 0.00038413 0.1343 999 0.928
## Residuals 105 0.300301 0.00286001#Homogeneity of Dispersion test with beta dispr
#Soerson DNA
beta_soren_station_DNA <- betadisper(scaled_sorensen_dist_DNA, metadata_DNA$lakesite)
permutest(beta_soren_station_DNA)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 0.06805 0.0226834 2.5481 999 0.056 .
## Residuals 52 0.46291 0.0089021
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Bray-curtis DNA
beta_bray_station_DNA <- betadisper(scaled_bray_dist_DNA, metadata_DNA$lakesite)
permutest(beta_bray_station_DNA)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 0.03891 0.0129706 1.6156 999 0.217
## Residuals 52 0.41748 0.0080285# Weighted Unifrac DNA
beta_bray_unifrac_DNA <- betadisper(scaled_wUnifrac_dist_DNA, metadata_DNA$lakesite)
permutest(beta_bray_unifrac_DNA)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 0.007897 0.0026323 1.1012 999 0.364
## Residuals 52 0.124307 0.0023905#Homogeneity of Dispersion test with beta dispr
#Soerson RNA
beta_soren_station_RNA <- betadisper(scaled_sorensen_dist_RNA, metadata_RNA$lakesite)
permutest(beta_soren_station_RNA)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 0.09634 0.032114 3.6571 999 0.02 *
## Residuals 49 0.43028 0.008781
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Bray-curtis RNA
beta_bray_station_RNA <- betadisper(scaled_bray_dist_RNA, metadata_RNA$lakesite)
permutest(beta_bray_station_RNA)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 0.07377 0.024588 2.3909 999 0.096 .
## Residuals 49 0.50393 0.010284
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Weighted Unifrac rNA
beta_bray_unifrac_RNA <- betadisper(scaled_wUnifrac_dist_RNA, metadata_RNA$lakesite)
permutest(beta_bray_unifrac_RNA)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 0.002733 0.00091112 0.2955 999 0.835
## Residuals 49 0.151078 0.00308323The high p-values from the betadisper test indicate that the variability is fairly homogeneous across our different samples.
Taxonomic Composition
#Phylum
# Set the phylum colors
phylum_colors <- c(
Acidobacteriota = "navy",
Actinobacteriota = "darkslategray2",
Armatimonadota = "deeppink1",
Alphaproteobacteria = "plum2",
Bacteroidota = "gold",
Betaproteobacteria = "plum1",
Bdellovibrionota = "red1",
Chloroflexi="black",
Crenarchaeota = "firebrick",
Cyanobacteria = "limegreen",
Deltaproteobacteria = "grey",
Desulfobacterota="magenta",
Firmicutes = "#3E9B96",
Gammaproteobacteria = "greenyellow",
"Marinimicrobia (SAR406 clade)" = "yellow",
Myxococcota = "#B5D6AA",
Nitrospirota = "palevioletred1",
Proteobacteria = "royalblue",
Planctomycetota = "darkorange",
"SAR324 clade(Marine group B)" = "olivedrab",
#Proteobacteria_unclassified = "greenyellow",
Thermoplasmatota = "green",
Verrucomicrobiota = "darkorchid1")
# Other = "grey")# Calculate the phylum relative abundance
# Make different dfs for combined, DNA and RNA samples
phylum_df_combined <-
scaled_rooted_physeq %>%
# agglomerate at the phylum level
tax_glom(taxrank = "Phylum") %>%
# Transform counts to relative abundance
transform_sample_counts(function (x) {x/sum(x)}) %>%
# Melt to a long format
psmelt() %>%
# Filter out Phyla that are <1% - get rid of low abundant Phyla
dplyr::filter(Abundance > 0.01)
phylum_df_DNA <-
scaled_rooted_physeq_DNA %>%
# agglomerate at the phylum level
tax_glom(taxrank = "Phylum") %>%
# Transform counts to relative abundance
transform_sample_counts(function (x) {x/sum(x)}) %>%
# Melt to a long format
psmelt() %>%
# Filter out Phyla that are <1% - get rid of low abundant Phyla
dplyr::filter(Abundance > 0.01)
phylum_df_RNA <-
scaled_rooted_physeq_RNA %>%
# agglomerate at the phylum level
tax_glom(taxrank = "Phylum") %>%
# Transform counts to relative abundance
transform_sample_counts(function (x) {x/sum(x)}) %>%
# Melt to a long format
psmelt() %>%
# Filter out Phyla that are <1% - get rid of low abundant Phyla
dplyr::filter(Abundance > 0.01)#Combined Phyla Plots
# Stacked Bar Plot With All phyla
# Plot Phylum Abundances - make sure to load phylum_colors
phylum_df_combined %>%
# Warning: It's important to have one sample per x value,
# otherwise, it will take the sum between multiple samples
dplyr::filter(limnion == "Bottom") %>%
ggplot(aes(x = lakesite, y = Abundance, fill = Phylum)) +
facet_grid(Phylum~season, scale = "free") +
geom_bar(stat = "identity", color = "black") +
labs(title = "Bottom Limnion Phylum Composition") +
scale_fill_manual(values = phylum_colors) +
theme_bw() +
theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1))
phylum_df_combined %>%
# Warning: It's important to have one sample per x value,
# otherwise, it will take the sum between multiple samples
dplyr::filter(limnion == "Middle") %>%
ggplot(aes(x = lakesite, y = Abundance, fill = Phylum)) +
facet_grid(Phylum~season, scale = "free") +
geom_bar(stat = "identity", color = "black") +
labs(title = "Middle Limnion Phylum Composition") +
scale_fill_manual(values = phylum_colors) +
theme_bw() +
theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1))
phylum_df_combined %>%
# Warning: It's important to have one sample per x value,
# otherwise, it will take the sum between multiple samples
dplyr::filter(limnion == "Top") %>%
ggplot(aes(x = lakesite, y = Abundance, fill = Phylum)) +
facet_grid(Phylum~season, scale = "free") +
geom_bar(stat = "identity", color = "black") +
labs(title = "Top Limnion Phylum Composition") +
scale_fill_manual(values = phylum_colors) +
theme_bw() +
theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1))
#Phyla plots from DNA samples
phylum_df_DNA %>%
# Warning: It's important to have one sample per x value,
# otherwise, it will take the sum between multiple samples
dplyr::filter(limnion == "Top") %>%
ggplot(aes(x = lakesite, y = Abundance, fill = Phylum)) +
facet_grid(Phylum~season, scale = "free") +
geom_bar(stat = "identity", color = "black") +
labs(title = "Top Limnion Phylum Composition DNA") +
scale_fill_manual(values = phylum_colors) +
theme_bw() +
theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1))
phylum_df_DNA %>%
# Warning: It's important to have one sample per x value,
# otherwise, it will take the sum between multiple samples
dplyr::filter(limnion == "Middle") %>%
ggplot(aes(x = lakesite, y = Abundance, fill = Phylum)) +
facet_grid(Phylum~season, scale = "free") +
geom_bar(stat = "identity", color = "black") +
labs(title = "Middle Limnion Phylum Composition DNA") +
scale_fill_manual(values = phylum_colors) +
theme_bw() +
theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1))
phylum_df_DNA %>%
# Warning: It's important to have one sample per x value,
# otherwise, it will take the sum between multiple samples
dplyr::filter(limnion == "Bottom") %>%
ggplot(aes(x = lakesite, y = Abundance, fill = Phylum)) +
facet_grid(Phylum~season, scale = "free") +
geom_bar(stat = "identity", color = "black") +
labs(title = "Bottom Limnion Phylum Composition DNA") +
scale_fill_manual(values = phylum_colors) +
theme_bw() +
theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1))
Top Limnion: For the total taxa based on lakesite we see slight increase in Actinobacteriota in summer throughout each lakesite and this seems to be the most abundant phyla. These are not the dominant active taxa, however. The dominant taxa according to the RNA samples are Verrucomicrobiota and Proteobacteria.Armatiomondota seem primarily to be at the surface and do not make appear in the active taxa. Cyanobacteria seem to be quite abundant in fall and summer across lake site. Interestingly the 15 meter lakesite seems to be unique in having cyanobateria in spring.
Middle Limnion: The only middle limnion lakesite sampled was 110m from the shore. The total taxa for the middle limnion is dominated by actinobacteriota and bacteroidota across season with marked increase in verrucomicrobiota in summer. This tracks well with the RNA samples, indicating similar trends between the total and active taxa.
Bottom Limnion: We similar trends of actinobacteriota to the top and middle limnion with similar trends across season. A point of note is that there seems to be no cyanobacteria present at any lakesite in spring. This limnion is also dominated by verrucomicriobiota. The active taxa track nicely with the total taxa for this depth, with the exception of Chloroflexi that make up a sizeable portion of the active taxa.
#Phyla from RNA samples
phylum_df_RNA <-
scaled_rooted_physeq_RNA %>%
# agglomerate at the phylum level
tax_glom(taxrank = "Phylum") %>%
# Transform counts to relative abundance
transform_sample_counts(function (x) {x/sum(x)}) %>%
# Melt to a long format
psmelt() %>%
# Filter out Phyla that are <1% - get rid of low abundant Phyla
dplyr::filter(Abundance > 0.01)
phylum_df_RNA %>%
# Warning: It's important to have one sample per x value,
# otherwise, it will take the sum between multiple samples
dplyr::filter(limnion == "Middle") %>%
ggplot(aes(x = lakesite, y = Abundance, fill = Phylum)) +
facet_grid(Phylum~season, scale = "free") +
geom_bar(stat = "identity", color = "black") +
labs(title = "Middle Limnion Phylum Composition RNA") +
scale_fill_manual(values = phylum_colors) +
theme_bw() +
theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1)) 
phylum_df_RNA %>%
# Warning: It's important to have one sample per x value,
# otherwise, it will take the sum between multiple samples
dplyr::filter(limnion == "Bottom") %>%
ggplot(aes(x = lakesite, y = Abundance, fill = Phylum)) +
facet_grid(Phylum~season, scale = "free") +
geom_bar(stat = "identity", color = "black") +
labs(title = "Bottom Limnion Phylum Composition RNA") +
scale_fill_manual(values = phylum_colors) +
theme_bw() +
theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1))
phylum_df_RNA %>%
# Warning: It's important to have one sample per x value,
# otherwise, it will take the sum between multiple samples
dplyr::filter(limnion == "Top") %>%
ggplot(aes(x = lakesite, y = Abundance, fill = Phylum)) +
facet_grid(Phylum~season, scale = "free") +
geom_bar(stat = "identity", color = "black") +
labs(title = "Top Limnion Phylum Composition RNA") +
scale_fill_manual(values = phylum_colors) +
theme_bw() +
theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1))
##Family
# Calculate the Family relative abundance for each sample
# # Make separate dfs for combined, DNA and RNA samples
family_df_combined <-
scaled_rooted_physeq %>%
# agglomerate at the phylum level
tax_glom(taxrank = "Family") %>%
# Transform counts to relative abundance
transform_sample_counts(function (x) {x/sum(x)}) %>%
# Melt to a long format
psmelt() %>%
# Filter out Phyla that are <1% - get rid of low abundant Phyla
dplyr::filter(Abundance > 0.01)
family_df_DNA <-
scaled_rooted_physeq_DNA %>%
# agglomerate at the phylum level
tax_glom(taxrank = "Family") %>%
# Transform counts to relative abundance
transform_sample_counts(function (x) {x/sum(x)}) %>%
# Melt to a long format
psmelt() %>%
# Filter out Phyla that are <1% - get rid of low abundant Phyla
dplyr::filter(Abundance > 0.01)
family_df_RNA <-
scaled_rooted_physeq_RNA %>%
# agglomerate at the phylum level
tax_glom(taxrank = "Family") %>%
# Transform counts to relative abundance
transform_sample_counts(function (x) {x/sum(x)}) %>%
# Melt to a long format
psmelt() %>%
# Filter out Phyla that are <1% - get rid of low abundant Phyla
dplyr::filter(Abundance > 0.01) #Judging from the Phylum plots we identified Cyanobacteria,Actinobacteriota and Verrucomicrobiota as the most interesting Phyla to discet further as they seem to be quite abudant and vary between season and limnion.
# Cyanobacteria Families in DNA sample
# Plot family
cyano_season_DNA <- family_df_DNA %>%
dplyr::filter(Phylum == "Cyanobacteria") %>%
# build the plot
ggplot(aes(x = lakesite, y = Abundance,
fill = lakesite, color = lakesite)) +
facet_grid(Family~season, scale = "free") +
geom_boxplot(alpha = 0.5, outlier.shape = NA) + # outliers not plotted here in boxplot
geom_jitter() +
theme_bw() +
labs(title = "Cyanobacteria Family Abundance") +
scale_color_manual(values = lakesite_colors) +
scale_fill_manual(values = lakesite_colors) +
theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1),
legend.position = "bottom")
actin_season_dna <- family_df_DNA %>%
dplyr::filter(Phylum == "Actinobacteriota") %>%
# build the plot
ggplot(aes(x = lakesite, y = Abundance,
fill = lakesite, color = lakesite)) +
facet_grid(Family~season, scale = "free") +
geom_boxplot(alpha = 0.5, outlier.shape = NA) + # outliers not plotted here in boxplot
geom_jitter() +
theme_bw() +
labs(title = "Actinobacteriota Family Abundance") +
scale_color_manual(values = lakesite_colors) +
scale_fill_manual(values = lakesite_colors) +
theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1),
legend.position = "bottom")
#actin_dna sanity check
verru_season_dna <- family_df_DNA %>%
dplyr::filter(Phylum == "Verrucomicrobiota") %>%
# build the plot
ggplot(aes(x = lakesite, y = Abundance,
fill = lakesite, color = lakesite)) +
facet_grid(Family~season, scale = "free") +
geom_boxplot(alpha = 0.5, outlier.shape = NA) + # outliers not plotted here in boxplot
geom_jitter() +
theme_bw() +
labs(title = " VerrucomicrobiotaFamily Abundance") +
scale_color_manual(values = lakesite_colors) +
scale_fill_manual(values = lakesite_colors) +
theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1),
legend.position = "bottom")
#verru_dna sanity check
#Now we will check the DNA families faceting by limnion
# Cyanobacteria Families in DNA sample
# Plot family
cyano_lim_DNA <- family_df_DNA %>%
dplyr::filter(Phylum == "Cyanobacteria") %>%
# build the plot
ggplot(aes(x = lakesite, y = Abundance,
fill = lakesite, color = lakesite)) +
facet_grid(Family~limnion, scale = "free") +
geom_boxplot(alpha = 0.5, outlier.shape = NA) + # outliers not plotted here in boxplot
geom_jitter() +
theme_bw() +
labs(title = "Cyanobacteria Family Abundance") +
scale_color_manual(values = lakesite_colors) +
scale_fill_manual(values = lakesite_colors) +
theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1),
legend.position = "bottom")
actin_lim_dna <- family_df_DNA %>%
dplyr::filter(Phylum == "Actinobacteriota") %>%
# build the plot
ggplot(aes(x = lakesite, y = Abundance,
fill = lakesite, color = lakesite)) +
facet_grid(Family~limnion, scale = "free") +
geom_boxplot(alpha = 0.5, outlier.shape = NA) + # outliers not plotted here in boxplot
geom_jitter() +
theme_bw() +
labs(title = "Actinobacteriota Family Abundance") +
scale_color_manual(values = lakesite_colors) +
scale_fill_manual(values = lakesite_colors) +
theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1),
legend.position = "bottom")
#actin_dna sanity check
verru_limnion_dna <- family_df_DNA %>%
dplyr::filter(Phylum == "Verrucomicrobiota") %>%
# build the plot
ggplot(aes(x = lakesite, y = Abundance,
fill = lakesite, color = lakesite)) +
facet_grid(Family~limnion, scale = "free") +
geom_boxplot(alpha = 0.5, outlier.shape = NA) + # outliers not plotted here in boxplot
geom_jitter() +
theme_bw() +
labs(title = " VerrucomicrobiotaFamily Abundance") +
scale_color_manual(values = lakesite_colors) +
scale_fill_manual(values = lakesite_colors) +
theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1),
legend.position = "bottom")
#verru_dna sanity check
#RNA Samples Plot
cyano_season_RNA <- family_df_RNA %>%
dplyr::filter(Phylum == "Cyanobacteria") %>%
# build the plot
ggplot(aes(x = lakesite, y = Abundance,
fill = lakesite, color = lakesite)) +
facet_grid(Family~season, scale = "free") +
geom_boxplot(alpha = 0.5, outlier.shape = NA) + # outliers not plotted here in boxplot
geom_jitter() +
theme_bw() +
labs(title = "Cyanobacteria RNA Family Abundance") +
scale_color_manual(values = lakesite_colors) +
scale_fill_manual(values = lakesite_colors) +
theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1),
legend.position = "bottom")
actin_season_rna <- family_df_RNA %>%
dplyr::filter(Phylum == "Actinobacteriota") %>%
# build the plot
ggplot(aes(x = lakesite, y = Abundance,
fill = lakesite, color = lakesite)) +
facet_grid(Family~season, scale = "free") +
geom_boxplot(alpha = 0.5, outlier.shape = NA) + # outliers not plotted here in boxplot
geom_jitter() +
theme_bw() +
labs(title = "Actinobacteriota Family RNA Abundance") +
scale_color_manual(values = lakesite_colors) +
scale_fill_manual(values = lakesite_colors) +
theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1),
legend.position = "bottom")
#actin_dna sanity check
verru_season_rna <- family_df_RNA %>%
dplyr::filter(Phylum == "Verrucomicrobiota") %>%
# build the plot
ggplot(aes(x = lakesite, y = Abundance,
fill = lakesite, color = lakesite)) +
facet_grid(Family~season, scale = "free") +
geom_boxplot(alpha = 0.5, outlier.shape = NA) + # outliers not plotted here in boxplot
geom_jitter() +
theme_bw() +
labs(title = " Verrucomicrobiota Family RNA Abundance") +
scale_color_manual(values = lakesite_colors) +
scale_fill_manual(values = lakesite_colors) +
theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1),
legend.position = "bottom")
#RNA Limnion Samples
cyano_lim_RNA <- family_df_RNA %>%
dplyr::filter(Phylum == "Cyanobacteria") %>%
# build the plot
ggplot(aes(x = lakesite, y = Abundance,
fill = lakesite, color = lakesite)) +
facet_grid(Family~limnion, scale = "free") +
geom_boxplot(alpha = 0.5, outlier.shape = NA) + # outliers not plotted here in boxplot
geom_jitter() +
theme_bw() +
labs(title = "Cyanobacteria RNA Family Abundance") +
scale_color_manual(values = lakesite_colors) +
scale_fill_manual(values = lakesite_colors) +
theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1),
legend.position = "bottom")
actin_lim_rna <- family_df_RNA %>%
dplyr::filter(Phylum == "Actinobacteriota") %>%
# build the plot
ggplot(aes(x = lakesite, y = Abundance,
fill = lakesite, color = lakesite)) +
facet_grid(Family~limnion, scale = "free") +
geom_boxplot(alpha = 0.5, outlier.shape = NA) + # outliers not plotted here in boxplot
geom_jitter() +
theme_bw() +
labs(title = "Actinobacteriota Family RNA Abundance") +
scale_color_manual(values = lakesite_colors) +
scale_fill_manual(values = lakesite_colors) +
theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1),
legend.position = "bottom")
#actin_dna sanity check
verru_lim_rna <- family_df_RNA %>%
dplyr::filter(Phylum == "Verrucomicrobiota") %>%
# build the plot
ggplot(aes(x = lakesite, y = Abundance,
fill = lakesite, color = lakesite)) +
facet_grid(Family~limnion, scale = "free") +
geom_boxplot(alpha = 0.5, outlier.shape = NA) + # outliers not plotted here in boxplot
geom_jitter() +
theme_bw() +
labs(title = " Verrucomicrobiota Family RNA Abundance") +
scale_color_manual(values = lakesite_colors) +
scale_fill_manual(values = lakesite_colors) +
theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1),
legend.position = "bottom")
cyano_lim_DNA
cyano_season_DNA
verru_limnion_dna
verru_season_rna
actin_season_dna
actin_season_rna
cyano_lim_RNA
verru_season_rna
actin_lim_rna
cyano_season_RNA
actin_season_rna
verru_season_rna
##Session Information
devtools::session_info()## ─ Session info ───────────────────────────────────────────────────────────────
## setting value
## version R version 4.3.2 (2023-10-31)
## os Rocky Linux 9.0 (Blue Onyx)
## system x86_64, linux-gnu
## ui X11
## language (EN)
## collate en_US.UTF-8
## ctype en_US.UTF-8
## tz America/New_York
## date 2024-04-30
## pandoc 3.1.1 @ /usr/lib/rstudio-server/bin/quarto/bin/tools/ (via rmarkdown)
##
## ─ Packages ───────────────────────────────────────────────────────────────────
## package * version date (UTC) lib source
## ade4 1.7-22 2023-02-06 [1] CRAN (R 4.3.2)
## ape 5.7-1 2023-03-13 [2] CRAN (R 4.3.2)
## Biobase 2.62.0 2023-10-24 [2] Bioconductor
## BiocGenerics 0.48.1 2023-11-01 [2] Bioconductor
## biomformat 1.30.0 2023-10-24 [1] Bioconductor
## Biostrings 2.70.1 2023-10-25 [2] Bioconductor
## bitops 1.0-7 2021-04-24 [2] CRAN (R 4.3.2)
## bslib 0.5.1 2023-08-11 [2] CRAN (R 4.3.2)
## cachem 1.0.8 2023-05-01 [2] CRAN (R 4.3.2)
## callr 3.7.3 2022-11-02 [2] CRAN (R 4.3.2)
## cli 3.6.1 2023-03-23 [2] CRAN (R 4.3.2)
## cluster 2.1.4 2022-08-22 [2] CRAN (R 4.3.2)
## codetools 0.2-19 2023-02-01 [2] CRAN (R 4.3.2)
## colorspace 2.1-0 2023-01-23 [2] CRAN (R 4.3.2)
## crayon 1.5.2 2022-09-29 [2] CRAN (R 4.3.2)
## data.table 1.14.8 2023-02-17 [2] CRAN (R 4.3.2)
## devtools * 2.4.4 2022-07-20 [2] CRAN (R 4.2.1)
## digest 0.6.33 2023-07-07 [2] CRAN (R 4.3.2)
## dplyr * 1.1.3 2023-09-03 [2] CRAN (R 4.3.2)
## ellipsis 0.3.2 2021-04-29 [2] CRAN (R 4.3.2)
## evaluate 0.23 2023-11-01 [2] CRAN (R 4.3.2)
## fansi 1.0.5 2023-10-08 [2] CRAN (R 4.3.2)
## farver 2.1.1 2022-07-06 [2] CRAN (R 4.3.2)
## fastmap 1.1.1 2023-02-24 [2] CRAN (R 4.3.2)
## forcats * 1.0.0 2023-01-29 [1] CRAN (R 4.3.2)
## foreach 1.5.2 2022-02-02 [2] CRAN (R 4.3.2)
## fs 1.6.3 2023-07-20 [2] CRAN (R 4.3.2)
## generics 0.1.3 2022-07-05 [2] CRAN (R 4.3.2)
## GenomeInfoDb 1.38.0 2023-10-24 [2] Bioconductor
## GenomeInfoDbData 1.2.11 2023-11-07 [2] Bioconductor
## ggplot2 * 3.5.0 2024-02-23 [2] CRAN (R 4.3.2)
## glue 1.6.2 2022-02-24 [2] CRAN (R 4.3.2)
## gridExtra * 2.3 2017-09-09 [1] CRAN (R 4.3.2)
## gtable 0.3.4 2023-08-21 [2] CRAN (R 4.3.2)
## highr 0.10 2022-12-22 [2] CRAN (R 4.3.2)
## hms 1.1.3 2023-03-21 [1] CRAN (R 4.3.2)
## htmltools 0.5.7 2023-11-03 [2] CRAN (R 4.3.2)
## htmlwidgets 1.6.2 2023-03-17 [2] CRAN (R 4.3.2)
## httpuv 1.6.12 2023-10-23 [2] CRAN (R 4.3.2)
## igraph 1.5.1 2023-08-10 [2] CRAN (R 4.3.2)
## IRanges 2.36.0 2023-10-24 [2] Bioconductor
## iterators 1.0.14 2022-02-05 [2] CRAN (R 4.3.2)
## jquerylib 0.1.4 2021-04-26 [2] CRAN (R 4.3.2)
## jsonlite 1.8.7 2023-06-29 [2] CRAN (R 4.3.2)
## knitr 1.45 2023-10-30 [2] CRAN (R 4.3.2)
## labeling 0.4.3 2023-08-29 [2] CRAN (R 4.3.2)
## later 1.3.1 2023-05-02 [2] CRAN (R 4.3.2)
## lattice * 0.21-9 2023-10-01 [2] CRAN (R 4.3.2)
## lifecycle 1.0.3 2022-10-07 [2] CRAN (R 4.3.2)
## lubridate * 1.9.3 2023-09-27 [1] CRAN (R 4.3.2)
## magrittr 2.0.3 2022-03-30 [2] CRAN (R 4.3.2)
## MASS 7.3-60 2023-05-04 [2] CRAN (R 4.3.2)
## Matrix 1.6-1.1 2023-09-18 [2] CRAN (R 4.3.2)
## memoise 2.0.1 2021-11-26 [2] CRAN (R 4.3.2)
## mgcv 1.9-0 2023-07-11 [2] CRAN (R 4.3.2)
## mime 0.12 2021-09-28 [2] CRAN (R 4.3.2)
## miniUI 0.1.1.1 2018-05-18 [2] CRAN (R 4.3.2)
## multtest 2.58.0 2023-10-24 [1] Bioconductor
## munsell 0.5.0 2018-06-12 [2] CRAN (R 4.3.2)
## nlme 3.1-163 2023-08-09 [2] CRAN (R 4.3.2)
## pacman 0.5.1 2019-03-11 [1] CRAN (R 4.3.2)
## patchwork * 1.2.0.9000 2024-04-28 [1] Github (thomasp85/patchwork@d943757)
## permute * 0.9-7 2022-01-27 [1] CRAN (R 4.3.2)
## phyloseq * 1.46.0 2023-10-24 [1] Bioconductor
## pillar 1.9.0 2023-03-22 [2] CRAN (R 4.3.2)
## pkgbuild 1.4.2 2023-06-26 [2] CRAN (R 4.3.2)
## pkgconfig 2.0.3 2019-09-22 [2] CRAN (R 4.3.2)
## pkgload 1.3.3 2023-09-22 [2] CRAN (R 4.3.2)
## plyr 1.8.9 2023-10-02 [2] CRAN (R 4.3.2)
## prettyunits 1.2.0 2023-09-24 [2] CRAN (R 4.3.2)
## processx 3.8.2 2023-06-30 [2] CRAN (R 4.3.2)
## profvis 0.3.8 2023-05-02 [2] CRAN (R 4.3.2)
## promises 1.2.1 2023-08-10 [2] CRAN (R 4.3.2)
## ps 1.7.5 2023-04-18 [2] CRAN (R 4.3.2)
## purrr * 1.0.2 2023-08-10 [2] CRAN (R 4.3.2)
## R6 2.5.1 2021-08-19 [2] CRAN (R 4.3.2)
## Rcpp 1.0.11 2023-07-06 [2] CRAN (R 4.3.2)
## RCurl 1.98-1.13 2023-11-02 [2] CRAN (R 4.3.2)
## readr * 2.1.5 2024-01-10 [1] CRAN (R 4.3.2)
## remotes 2.4.2.1 2023-07-18 [2] CRAN (R 4.3.2)
## reshape2 1.4.4 2020-04-09 [2] CRAN (R 4.3.2)
## rhdf5 2.46.1 2023-11-29 [1] Bioconductor 3.18 (R 4.3.2)
## rhdf5filters 1.14.1 2023-11-06 [1] Bioconductor
## Rhdf5lib 1.24.2 2024-02-07 [1] Bioconductor 3.18 (R 4.3.2)
## rlang 1.1.2 2023-11-04 [2] CRAN (R 4.3.2)
## rmarkdown 2.25 2023-09-18 [2] CRAN (R 4.3.2)
## rstudioapi 0.15.0 2023-07-07 [2] CRAN (R 4.3.2)
## S4Vectors 0.40.1 2023-10-26 [2] Bioconductor
## sass 0.4.7 2023-07-15 [2] CRAN (R 4.3.2)
## scales 1.3.0 2023-11-28 [2] CRAN (R 4.3.2)
## sessioninfo 1.2.2 2021-12-06 [2] CRAN (R 4.3.2)
## shiny 1.7.5.1 2023-10-14 [2] CRAN (R 4.3.2)
## stringi 1.7.12 2023-01-11 [2] CRAN (R 4.3.2)
## stringr * 1.5.0 2022-12-02 [2] CRAN (R 4.3.2)
## survival 3.5-7 2023-08-14 [2] CRAN (R 4.3.2)
## tibble * 3.2.1 2023-03-20 [2] CRAN (R 4.3.2)
## tidyr * 1.3.0 2023-01-24 [2] CRAN (R 4.3.2)
## tidyselect 1.2.0 2022-10-10 [2] CRAN (R 4.3.2)
## tidyverse * 2.0.0 2023-02-22 [1] CRAN (R 4.3.2)
## timechange 0.3.0 2024-01-18 [1] CRAN (R 4.3.2)
## tzdb 0.4.0 2023-05-12 [1] CRAN (R 4.3.2)
## urlchecker 1.0.1 2021-11-30 [2] CRAN (R 4.3.2)
## usethis * 2.2.2 2023-07-06 [2] CRAN (R 4.3.2)
## utf8 1.2.4 2023-10-22 [2] CRAN (R 4.3.2)
## vctrs 0.6.4 2023-10-12 [2] CRAN (R 4.3.2)
## vegan * 2.6-4 2022-10-11 [1] CRAN (R 4.3.2)
## withr 2.5.2 2023-10-30 [2] CRAN (R 4.3.2)
## xfun 0.41 2023-11-01 [2] CRAN (R 4.3.2)
## xtable 1.8-4 2019-04-21 [2] CRAN (R 4.3.2)
## XVector 0.42.0 2023-10-24 [2] Bioconductor
## yaml 2.3.7 2023-01-23 [2] CRAN (R 4.3.2)
## zlibbioc 1.48.0 2023-10-24 [2] Bioconductor
##
## [1] /home/dt473/R/x86_64-pc-linux-gnu-library/4.3
## [2] /programs/R-4.3.2/library
##
## ──────────────────────────────────────────────────────────────────────────────